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ABSTRACT 

We propose a scenario that explains the apparent nitrogen deficiency in 
comets in a way consistent with the fact that the surfaces of Pluto and 
Triton are dominated by nitrogen-rich ice. We use a statistical thermo- 
dynamic model to investigate the composition of the successive multiple guest 
clathrates that may have formed during the cooling of the primordial nebula 
from the most abundant volatiles present in the gas phase. These clathrates 
agglomerated with the other ices (pure condensates or stoichiometric hydrates) 
and formed the building blocks of comets. We report that molecular nitrogen is 
a poor clathrate former, when we consider a plausible gas phase composition of 
the primordial nebula. This implies that its trapping into cometesimals requires 
a low disk temperature (~20 K) in order to allow the formation of its pure con- 
densate. We find that it is possible to explain the lack of molecular nitrogen in 
comets as a consequence of their postformation internal heating engendered by 
the decay of short-lived radiogenic nuclides. This scenario is found consis- 
tent with the presence of nitrogen-rich ice covers on Pluto and Triton. Our 
model predicts that comets should present xenon-to-water and krypton-to-water 
ratios close to solar xenon-to-oxygen and krypton-to-oxygen ratios, respectively. 
In contrast, the argon-to- water ratio is predicted to be depleted by a factor of 
~300 in comets compared to solar argon-to-oxygen, as a consequence of poor 
trapping efficiency and radiogenic heating. 

Subject headings: Comets: general - Kuiper Belt: general - Protoplanetary disks - 
Astrobiology 
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1. Introduction 

The abundance of the super volatile molecule N2 in cometary nuclei is still under 
debate. Searches for UV fluorescence from N 2 with the Far Ultraviolet Spectroscopic 
Explorer (FUSE) have been unsuccessful (Bockelee-Morvan et al. 2004). Measurements of 
the N2 abundance in lP/Halley by Giotto were not feasible: owing to the low resolution of 
the mass spectrometer, it was not possible to discriminate N2 from CO (Eberhardt et al. 
1987). Ground-based observation of the band at 3914 A is a difficult challenge due both 
to the contamination of the C3 (0,2,0)-(0,0,0) band and the presence of telluric emission 
lines. Different positive detections of this feature have been claimed in the ionic tail of 
comets Humason (1962 VIII) (Greenstein & Arpigny 1962), lP/Halley (Lutz et al. 1993; 
Wyckoff & Theobald 1989), C/1987 PI Bradfield (Lutz et al. 1993), and C/2002 VQ94 
(LINEAR) (Korsun et al. 2006). Arpigny, from photographic spectra and other archival 
data estimated a positive N^/CO + intensity ratio for 12 comets (Cochran et al. 2000). 
All these positive detections are based on low-resolution spectra. So far similar searches 
with high-resolution spectra have been unsuccessful. This was the case for observations of 
comets 122P/de Vico, C/1995 01 (Hale-Bopp), and 153P/2002 CI (Ikeya-Zhang) that did 
not detect band, yielding upper limits of < 10~ 5 -10~ 4 on the abundance of N2 relative 
to CO (Cochran et al. 2000; Cochran 2002). 

An interpretation of the possible N2 deficiency in comets has been proposed by Iro et 
al. (2003). These authors assumed that comets were made of clathrate hydrates (hereafter 
clathrates) and argued that CO was preferentially incorporated in clathrates compared to 
N 2 after having investigated the competition between the trapping of these two molecules 
in clathrates. In order to explain the N2 deficiency in comets, they argued that the nebula's 
temperature never cooled down below ~45 K in their formation region, impeding the 
formation of N 2 -bearing ice (in clathrate or pure condensate form) that would have been 
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agglomerated by comets. However, this mechanism is not consistent with the fact that Pluto 
and Triton possess thick nitrogen ice covers (Lellouch et al. 2011) while they are both 
expected to have been formed in the same region of the primitive nebula as ecliptic comets 
(Kavelaars et al. 2001). In addition, Iro et al. (2003) did not consider the competition 
between the trapping of N2 in clathrates with that of the other abundant volatiles present in 
the nebula (CO, C0 2 , CH 4 , H 2 S, ...). Indeed, due to their various propensities for trapping, 
the consideration of a larger set of molecules can drastically change the calculated fraction 
of N 2 incorporated in clathrates formed in the nebula. 

Here we propose a scenario that could explain the apparent N 2 deficiency in comets 
in a way consistent with the presence of this molecule on the surfaces of Pluto and 
Triton. We use a statistical thermodynamic model to investigate the composition of the 
successive multiple guest clathrates (hereafter MG clathrates) that may have formed during 
the cooling of the primordial nebula from the most abundant volatiles present in the gas 
phase (see section 2). These clathrates agglomerated with the other ices and formed the 
building blocks of comets (hereafter cometesimals) . The major ingredient of our model is 
the description of the guest-clathrate interaction by a spherically averaged Kihara potential 
with a nominal set of potential parameters. Our model allows us to report that N 2 is a 
poor clathrate former when considering a plausible gas phase composition of the primordial 
nebula, implying that its trapping into cometesimals requires a low disk's temperature 
(~20 K) in order to allow the formation of its pure condensate. We find that it is possible 
to explain the lack of nitrogen in comets as a consequence of their postaccretion internal 
heating engendered by the decay of radiogenic nuclides (see section 3). 
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2. Composition of cometesimals formed in the primordial nebula 

2.1. Modeling approach 

In our model, the volatile phase incorporated in cometesimals is composed of a 
mixture of pure ices, stoichiometric hydrates (such as NH 3 -H 2 hydrate) and MG 
clathrates that crystallized in the form of microscopic grains at various temperatures in 
the outer part of the disk. Our model is based on the assumption that cometesimals have 
grown from the agglomeration of these icy grains due to collisional coagulation, implying 
no loss of volatile during their growth phase (Weidenschilling 1997). Here, the clathration 
process stops when no more crystalline water ice is available to trap the volatile species 
and then only pure condensates form at lower temperature in the disk. The process of 
volatiles trapping in icy grains is calculated using the equilibrium curves of hydrates and 
pure condensates, our model determining the equilibrium curves and compositions of 
MG clathrates, and the thermodynamic paths detailing the evolution of temperature and 
pressure between 5 and 30 AU in the protoplanetary disk. We refer to the work of Alibert 
et al. (2005) for a full description of the turbulent model of accretion disk used here. 

The composition of the initial gas phase of the disk is defined as follows: we assume 
that the abundances of all elements (C, N, O, S, P, Ar, Kr and Xe) are protosolar (Asplund 
et al. 2009) and that O, C, and N exist only under the form of H 2 0, CO, C0 2 , CH 3 OH, 
CH 4 , N 2 , and NH 3 . The abundances of CO, C0 2 , CH 3 OH, CH 4 , N 2 and NH 3 are then 
determined from the adopted CO:C0 2 :CH 3 OH:CH4 and N 2 :NH 3 gas phase molecular ratios. 
Once the abundances of these molecules are fixed, the remaining O gives the abundance 
of H 2 0. Concerning the distribution of elements in the main volatile molecules, we set 
CO:C0 2 :CH 3 OH:CH4 = 70:10:2:1 in the gas phase of the disk, values that are consistent 
with the ISM measurements considering the contributions of both gas and solid phases in 
the lines of sight (Mousis et al. 2009). In addition, S is assumed to exist in the form of H 2 S, 
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with H2S:H2 = 0.5 x(S/H2)©, and other refractory sulfide components. We also consider 
N 2 :NH 3 = 1:1 in the nebula gas-phase, a value compatible with thermochemical calculations 
in the solar nebula that take into account catalytic effects of Fe grains on the kinetics of N 2 
to NH3 conversion (Fegley 2000). Note that CH3OH is considered to be formed only as a 
pure condensate in our calculations since, to the best of our knowledge, no experimental 
data concerning the equilibrium curve of its associated clathrate have been reported in the 
literature. 

The equilibrium pressure for the successive MG clathrates formed in the solar nebula 
can be expressed as (Hand et al. 2006): 



where yi is the mole fraction of the component i in the fluid phase. The equilibrium pressure 
curves of each species are determined by fitting the available theoretical and laboratory 
data (Lunine & Stevenson 1985). Their equations are of the form \ogP eq ^ = A/T + B, 
where P eq ^ and T are the partial equilibrium pressure (bars) and temperature (K) of the 
considered species i, respectively. 

The relative abundances of guest species incorporated in MG clathrate formed at a 
given temperature and pressure from the solar nebula gas phase are calculated following 
the method described by Lunine & Stevenson (1985) and Mousis et al. (2010), which 
uses classical statistical mechanics to relate the macroscopic thermodynamic properties of 
clathrates to the molecular structure and interaction energies. It is based on the original 
ideas of van der Waals & Platteeuw (1959) for clathrate formation, which assume that 
trapping of guest molecules into cages corresponds to the three-dimensional generalization 
of ideal localized adsorption. This approach is based on four key assumptions (Lunine & 
Stevenson 1985; Sloan 1998): 
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1. The host molecules contribution to the free energy is independent of the clathrate 
occupancy. This assumption implies in particular that the guest species do not distort 
the cages. 

2. (a) The cages are singly occupied, (b) Guest molecules rotate freely within the cage. 

3. Guest molecules do not interact with each other. 

4. Classical statistics is valid, i.e., quantum effects are negligible. 

In this formalism, the fractional occupancy of a guest molecule K for a given type t 
(t = small or large) of cage can be written as 



where the sum in the denominator includes all the species which are present in the initial 
gas phase. Cx,t is the Langmuir constant of species K in the cage of type t, and Pk is the 
partial pressure of species K. This partial pressure is given by Pk = xk x P (we assume 
that the sample behaves as an ideal gas), with xk the mole fraction of species K in the gas 
phase, and P the total gas pressure, which is dominated by H 2 . 

The Langmuir constant depends on the strength of the interaction between each guest 
species and each type of cage, and can be determined by integrating the molecular potential 
within the cavity as 



where R c represents the radius of the cavity assumed to be spherical, ks the Boltzmann 
constant, and WK,t{ r ) is the spherically averaged Kihara potential representing the 



VK,t = 




(2) 
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interactions between the guest molecules K and the H2O molecules forming the surrounding 
cage t. For a spherical guest molecule, this potential w(r) can be written as (McKoy & 
Sinanoglu 1963): 
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with 



In Eq. (jlj), z is the coordination number of the cell. This parameter depends on the 
structure of the clathrate (I or II) and on the type of the cage (small or large). The Kihara 
parameters a, a and e for the molecule- water interactions employed in this work are listed 
in Table [TJ Note that our parameters are different from those used by Iro et al. (2003), 
who instead of using potential parameters describing the guest-clathrate interaction in their 
statistical model, employed potential parameters corresponding to guest-guest interactions 
(case of pure solutions). When comparing different data sets determined for the same 
X-H2O interaction (with X = CH4, N2 and Xe), we opted for the most recent ones because 
they are fitted to a larger range of experimental data. In the case of Ar, Kr, CO and 
PH 3 , the listed sets are the unique ones that we found in the literature. For the H 2 S-H 2 
and CO2-H2O interaction parameters we used those of Parrish and Prausnitz (1972) and 
Thomas et al. (2009) (the set of CO2-H2O interaction parameters provided by Thomas et 
al. (2009) is an update of the one determined by Parrish and Prausnitz (1972) and very 
close to the one recently derived by Herri and Chassefiere (2012)). We also investigated 
the interaction parameters of Sloan and Koh (2008) but these parameters yield unphysical 
results when used in our models. We suspect this is because the code of Sloan and Koh 
(2008) is optimized for industrial purposes and the parameters are not appropriate for 
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our temperature-pressure regime. We have tested our code against the one of Herri and 
Chassefiere (2012) so we do not believe the problem resides in our code. Finally, the mole 
fraction fx of a guest molecule K in a clathrate can be calculated with respect to the whole 
set of species considered in the system as 



, b s y K , s + b e y K ,e 

where b s and bi are the number of small and large cages per unit cell respectively, for 
the clathrate structure under consideration, and with Jk — 1. Values of z, R c , b s 
and bi are taken from Parrish & Prausnitz (1972). Each time a species is fully trapped 
in a MG clathrate, its gaseous abundance is set to zero in the disk. This change of the 
disk's gas phase composition induces the formation of a new MG clathrate at different 
temperature and pressure conditions and with a distinct composition. As a consequence, 
the compositions and formation conditions of successive MG clathrates are calculated n 
times in the nebula, with n corresponding to the total number of fully enclathrated species. 
Once all the water budget has been used for clathration, the volatiles remaining in the gas 
phase form pure condensates at lower temperatures. The equilibrium curves of these pure 
condensates derive from the compilation of laboratory data given in the CRC Handbook of 
Chemistry and Physics (Lide 2002). 

The composition of the volatile phase incorporated into cometesimals formed at a 
given distance from the Sun is finally given by the intersection of the disk's thermodynamic 
path at this location with the equilibrium curves of the different formed ices. The amount 
of volatile, i (relative to water), used to form the ice j in the nebula (either in the form 
of stoichiometric hydrate, MG clathrate or of pure condensate), can be determined by the 
relation: 
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where is the mass mixing ratio with respect to H2 of the volatile i used to form ice j 
in the nebula. Xh 2 o is the mass mixing ratio of H2O with respect to H2 in the nebula. 
Tl(R;Tj, Pj) and Y<(R]T H2 o, Ph 2 o) are the surface density of the disk at distance R from 
the Sun at the epoch of formation of ice j, and at the epoch of condensation of water, 
respectively. The global volatile, i, to water mass ratio incorporated in cometesimals is 
then given by Mj = 52/=i k m v with k corresponding to the number of solid phases in 
which volatile % is incorporated. Note that a thermodynamic path has been arbitrarily 
selected at the heliocentric distance of 5 AU to determine the composition of the formed 
ices. The adoption of any other distance range between ~5 and 30 AU in the nebula would 
not affect the composition of the ices because it remains almost identical irrespective of i) 
their formation distance and ii) the input parameters of the disk, provided that the initial 
gas phase composition is homogeneous (Marboeuf et al. 2008). 



Figure [T] represents the composition of the volatile phase incorporated in cometesimals 
agglomerated from ices condensed in the outer nebula and Figure [2] displays the molar 
ratios between incorporated volatiles of interest. Both figures show calculations expressed 
as a function of the formation temperature of cometesimals in the primordial nebula. 
Observations of active comets when their activity is driven by the sublimation of water 
ice - i.e. for heliocentric distances smaller than about 2 AU - suggest that the mean CO 
production rate relative to water is typically a few percents, with extreme values comprised 
between 0.4 and 30% (Bockelee-Morvan et al. 2004). Assuming that the production 
rates relative to water measured in comets are representative of the bulk composition of 
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Results 
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cometesimals, Figure [2] shows that their formation temperature must be lower than ~47 
K in the nebula to get CO/H 2 > 1% in their interior. The figure also shows that a 
formation temperature lower than ~25K would increase the CO/H 2 ratio up to ~13% in 
cometesimals, as a result of the incorporation of pure CO condensate at this lower nebular 
temperature. 

Formation temperatures of ices lower than ~22 K in the nebula allow trapping of N 2 
with N2/H2O of order of 1% in the interior of planetesimals, a value far above the highest 
estimates (~0.02%) inferred from observations (Wyckoff et al. 1991; Cochran et al. 2000). 
On the other hand, assuming a formation temperature of cometesimals in the 22-47 K 
range in the nebula would allow N2/H2O ~0.01% and CO/H2O ~6.6% in their interior, 
values that are consistent with observations. A simple explanation for the formation of 
cometesimals would then argue that they need to be agglomerated from grains formed in 
the 22-47 K range in order to match the observations of the cometary production rates. 
However, this mechanism is not consistent with the fact that Pluto and Triton possess 
N2-rich ice covers while they are expected to be formed in the same region as comets. 
An alternative possibility would be that comets accreted from planetesimals formed at very 
low temperature in the disk, but that subsequent radiogenic heating enabled their partial 
devolatilization through their pores. In the next Section, we explore how radiogenic heating 
can induce important losses of N 2 in comets after their formation. 



3. Thermal evolution of comets 
3.1. Nucleus model 



We use a fully thr e e dim ensional model of heat transport described by 
Guilbert-Lepoutre et al.l (120111 ) to compute the time evolution of the temperature 
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distribution within cometary nuclei, accounting for the heating due to the radioactive 
decay of short-lived nuclides 26 Al and 60 Fe. We assume that these objects formed by cold 
accretion, i.e. no accretional heating is accounted for. The model computes the temperature 
distribution and its evolution with time, by solving the heat diffusion equation: 



PbulkC 



dT 
~dt 



V(-k^T) = Q 



rad 



with pbuik [kg m -3 ] the bulk density, c [J kg -1 K _1 ] the heat capacity, k [W hi -1 K _1 ] the 
thermal conductivity and Q ra d [W m -3 ] the internal power production per unit volume due 
to the decay of short-lived radioactive nuclides. As an exa mple, we show in this paper the 



resul t s for an idealized "Hale-Bopp" wi ; 
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2000 



Weaver fc Lamy 



i a s ize of 35 km fjSzabo et al. 



2011 
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19971 ). We consider this object sphere, made 



of a porous mixture of crystalline water ice and dust, homogeneously distributed within 
the ice matrix. We use generic parameters to describe the bulk material: we consider it 
can be described with a dust to water ice mass ratio Xd/X w = 1 (with and X w the 
mass fractions of dust and water respectively), a porosity ip = 50%, and a bulk de nsity 



Huebner et al 



(|2006|)). The 



Pbuik = 700 kg m~ 3 (typical values of the comet literature; see 
heat capacity and thermal conductivity of such a mixture are c = 800 J kg -1 K _1 and k = 
0.09 W m -1 K" 1 respectively. The thermal evolution is computed with the formation time 
tp -time for an object to grow to its final size- as a reference for time zero, from an initial 
temperature of 20 K. 



Heating due to radioactive decay is described by: 

Qrad 



^Pd X rad (t F ) H rad —exp ( 

rad Trad V 



-t 

Trad 



(9) 



with pd the dust bulk density, X ra( i{tp) the mass fraction of the given radioactive isotope after 
a formation time tp, H rad the heat released per unit mass upon decay (4.84 x 10 12 J kg -1 
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and 5.04 x 10 12 J kg 1 for 26 Al and 60 Fe respectively (see iGuilbert-Lepoutre et al. 



T ra d its mean-lifetime (using 1.05 Myr for 26 Al and 3.78 Myr for 60 Fe f lRugel et al. 



20111) 



2009) 



and t the time. The object formation time tp is accounted for in the model through the 
decay of radioactive isotopes that occurs during accretion. No additional accretional heating 
is accounted for, though. Decay during the formation time t F results in a decrease of each 
nuclide initial abundance in the simulations: 



X rad {t F ) = X rad (0) e 



-tp/T ra d 



(10) 



with X ra d(0) the initial mass fraction of each nuclide. We consider the initial nuclide 



ratios 26 A1/ 27 A1 = 5 x 10" 5 ( IMacPherson et al 



( iSpivak-Birndorf et al. 



20111 ). and mass fractions from 



19951) and 60 Fe/ 56 Fe = 1 x 10 



Wasson & KallemevrJ (119881 ) for CV 



chondrites. 



3.2. Sensivity to parameters 

Although the model has many free parameters, only a few can actually strongly affect 
the resulting thermal history, such as the formation time of the object, its composition 
or the thermal conductivity. We also point out that the object's size plays an important 
role in its thermal evolution, owing to the surface to volume ratio. Figure [3] represents the 
influence of the formation time, through the abundance of radiogenic nuclides. It shows 
the evolution of the object's central temperature as a function of time after formation, 
for different formation times. The central temperature corresponds to the maximum 
temperature achieved inside the object: since the surface is in thermal equilibrium with 
the surrounding nebula at these stages of evolution, the heating due to radioactive decay 
is balanced by a cold load coming from the surface. Figure H] displays the time evolution 
of the temperature profile within the body, calculated from a formation delay of 3 Myr. 
This figure shows that, at a given epoch of its evolution and for the adopted set of input 
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parameters, the condition that the temperature profile must be in the 22-47 K range can be 
satisfied in almost the whole body, thus enabling the release of N 2 from this zone. Only the 
very close subsurface (down to a few hundreds meters deep) remains at lower temperatures, 
but this region should be rapidly depleted in volatiles once the comet follows subsequent 
orbits closer to the Sun. 

The combined effects of composition and porosity are presented in Figure [5j For those 
simulations, we also considered a plausible formation time of 3 Myr and a fixed density 
of 700 kg m -3 . Figure M represents the influence of thermal conductivity for the standard 
parameters previously mentioned. All together, these simulations show that it is possible to 
heat planetesimals at temperatures high enough to simultaneously imply the loss of N 2 and 
the preservation of CO, between 22 and 47 K. 

4. Discussion 

Our computations show that, under certain sets of conditions and 
evolutionary paths, initial nitrogen-rich cometesimals similar to Triton and 
Pluto may have evolved into nitrogen-depleted comets. The proposed scenario 
of selective trapping in clathrates in the nebula followed by outgassing due to 
short-lived radiogenic heating thus reconciles present-day observations with 
solar system formation models. In our scenario, the composition of comets would have 
evolved as a result of radiogenic decay of short-lived nuclides. In the case of large objects 
like Pluto or Triton, the thermal evolution is more complex and involves partial or total 
physical differentiation. Although this significant processing could involve a substantial loss 
of volatiles, Pluto and Triton are actually large enough for most volatiles to be retained 
at the surface and in the atmosphere by gravity, as shown by Schaller & Brown (2007). 
Moreover, if accretion of these bodies were occurring in a cold environment and were 
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slow-so that radiative removal of the impact heat were effective- then Triton, Pluto and 
other large KBOs could have retained their volatiles and accreted fairly cold (see discussion 
by Canup & Ward (2002) in the context of the accretion of Ganymede and Callisto). 

Our model predicts that comets should contain Xe/H^O (~3.6 x 10~ 7 ) and Kr/B^O 
(~5.0 x 10~ 6 ), ratios close to solar Xe/O and Kr/O, respectively. In contrast, the Ar/H 2 
ratio is predicted to be equal to ~1.7 x 10~ 5 in comets, a value about 300 times lower than 
solar Ar/O. In our model, the Ar depletion (compared to solar) in comets results from the 
poor propensity of this element to be trapped in different MG clathrates formed in the 
nebula. Because of this property, Ar was accreted by cometesimals only at low nebular 
temperature (~20 K) and in pure condensate form. As a consequence, the radiogenic 
heating above 22 K triggers the dissociation of clathrates and the outgassing of 
N2 and sublimation of Ar. 

Measuring the noble gas abundances in comets might require in situ measurements via 
a high sensitivity mass spectrometer or ultraviolet spectroscopy. The ALICE ultraviolet 
spectrometer and the ROSINA instrument should put upper limits on the noble gas 
abundances in the coma of comet 67P/Churyumov-Gerasimenko (Balsiger et al. 2007; 
Stern et al. 2007) that will probably help understand the formation conditions of 
comets in the primordial nebula. However, the measurement of solar or subsolar noble gas 
abundances in comets will require a new generation of instruments such as the cryotrap, 
which is part of the mass spectrometer MASPEX developed in the context of the PRIME 
spacecraft project (Young et al. 2010). 

O. M. and J.-M. P. acknowledge support from CNES, and J. L. from the Distinguished 
Visiting Scientist Program at JPL. We acknowledge an anonymous Referee for his 
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Table 1: Adopted set of parameters for the Kihara potential 



Mol. 


a (A.) 


e/k B (K) 


a (A) 


Reference 


CO 


3.1515 


133.61 


0.3976 


Mohammadi et al. (2005) 


co 2 


2.9681 


169.09 


0.6805 


Thomas et al. (2009) 


CH 4 


3.14393 


155.593 


0.3834 


Sloan & Koh (2008) 


N 2 


3.13512 


127.426 


0.3526 


Sloan & Koh (2008) 


H 2 S 


3.1558 


205.85 


0.36 


Parrish & Prausnitz (1972) 


PH 3 


3.771 


275.0 





Vorotyntsev & Malyshev (1998) 


Ar 


2.9434 


170.50 


0.184 


Parrish & Prausnitz (1972) 


Kr 


2.9739 


198.34 


0.230 


Parrish & Prausnitz (1972) 


Xe 


3.32968 


193.708 


0.2357 


Sloan & Koh (2008) 



a is the Lennard- Jones diameter, e is the depth of the potential well, and a is the radius of 
the impenetrable core. 
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Fig. 1. — Composition of the volatile phase incorporated in cometesimals as a function of 
their formation temperature in the outer solar nebula. 
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Fig. 2. — Molar ratios between different volatiles calculated in cometesimals as a function of 
their formation temperature. 
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Fig. 3. — Time evolution of central temperature calculated fror different formation times. 
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Fig. 4. — Time evolution of the temperature profile within the body after a formation delay 
of 3 Myr. The parameters of the model are those adopted in Fig. [31 Top panel represents 
a global view of the body's heating due to radionuclides and then cooling over a timespan 
of 85 Myr after formation. Bottom panel represents a zoom of the temperature evolution of 
the body's subsurface over 2.5 Myr after formation. 
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Fig. 5. — Time evolution of central temperature calculated for different values of the porosity 
and composition. 
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Fig. 6. — Time evolution of central temperature after a formation delay of 3 Myr, calculated 
for different values of the thermal conductivity. 



